17.2: Robust regression

Introduction

This is code for Section 17.2 of Kruschke using a Student t noise distribution. Other than that change, this document is very similar to the one from section 7.1.

Make sure that the two files “HtWtData30.csv” and “HtWtData300.csv” are in your project directory.

Preliminaries

Load necessary packages:

library(tidyverse)
library(rstan)
library(tidybayes)
library(bayesplot)

Set Stan to save compiled code.

rstan_options(auto_write = TRUE)

Set Stan to use parallel processing where possible.

options(mc.cores = parallel::detectCores())

Data

HtWtData30 <- read_csv("HtWtData30.csv")

Rows: 30 Columns: 3── Column specification ─────────────────────────────────────────
Delimiter: ","
dbl (3): male, height, weight
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

HtWtData30

ggplot(HtWtData30, aes(y = weight, x = height)) +
    geom_point()

HtWtData300 <- read_csv("HtWtData300.csv")

Rows: 300 Columns: 3── Column specification ─────────────────────────────────────────
Delimiter: ","
dbl (3): male, height, weight
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

HtWtData300

ggplot(HtWtData300, aes(y = weight, x = height)) +
    geom_point()

Mean center the x variable. (Kruschke standardizes both variables using z-scores, but that makes the coefficients of the model harder to interpret.)

HtWtData30 <- HtWtData30 %>%
    mutate(height_mc = height - mean(height))
HtWtData30

ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
    geom_point()

HtWtData300 <- HtWtData300 %>%
    mutate(height_mc = height - mean(height))
HtWtData300

ggplot(HtWtData300, aes(y = weight, x = height_mc)) +
    geom_point()

We’re going to run two different examples, so we’ll make two lists.

stan_data_30 <- list(N = NROW(HtWtData30),
                     y = HtWtData30$weight,
                     x = HtWtData30$height_mc)
stan_data_300 <- list(N = NROW(HtWtData300),
                      y = HtWtData300$weight,
                      x = HtWtData300$height_mc)

Prior predictive check

Simulation code

data {
    int<lower = 0> N;
    vector<lower = 0>[N] y;
    vector[N] x;
}
generated quantities {
    real nu;
    real beta0;
    real beta1;
    vector[N] mu;
    real<lower = 0> sigma;
    real y_sim[N];
    
    nu = gamma_rng(2, 0.1);       // default prior for df
    beta0 = normal_rng(150, 25);  // Intercept between 100 and 200
    beta1 = normal_rng(0, 5);    // -10 to 10 lbs per inch
    mu = beta0 + beta1 * x;       // linear model
    sigma = exponential_rng(0.1); // default prior for sigma
    y_sim = student_t_rng(nu, mu, sigma);
}

Simulate data

fit_HtWt30_robust_prior <- sampling(HtWt_robust_prior,
                                    data = stan_data_30,
                                    chains = 1,
                                    algorithm = "Fixed_param",
                                    seed = 11111,
                                    refresh = 0)

Extract samples

samples_HtWt30_prior <- tidy_draws(fit_HtWt30_robust_prior)
samples_HtWt30_prior

y_sim_prior <- samples_HtWt30_prior %>%
    select(starts_with("y_sim")) %>%
    as.matrix()

Plot results

ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
    geom_point() +
    geom_abline(data = samples_HtWt30_prior,
                aes(intercept = beta0, slope = beta1),
                alpha = 0.05)

mcmc_hist(samples_HtWt30_prior,
          pars = c("nu", "sigma", "beta0", "beta1"))

ppc_intervals(y = HtWtData30$weight,
              x = HtWtData30$height,
              yrep = y_sim_prior)

Model

data {
    int<lower = 0> N;
    vector<lower = 0>[N] y;
    vector[N] x;
}
parameters {
    real nu;
    real beta0;
    real beta1;
    real<lower = 0> sigma;
}
transformed parameters {
    vector[N] mu;
    mu = beta0 + beta1 * x;   // linear model
}
model {
    nu ~ gamma(2, 0.1);       // default prior for df
    beta0 ~ normal(150, 25);  // Intercept between 100 and 200
    beta1 ~ normal(0, 5);    // -10 to 10 lbs per inch
    sigma ~ exponential(0.1); // default prior for sigma
    y ~ student_t(nu, mu, sigma);
}
generated quantities {
    real y_sim[N];
    y_sim = student_t_rng(nu, mu, sigma);
}

fit_HtWt30_robust <- sampling(HtWt_robust,
                              data = stan_data_30,
                              seed = 11111,
                              refresh = 0)

Warning: There were 676 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problems
Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

Diagnosing the model

plot(fit_HtWt30_robust, plotfun = "ac",
     pars = c("nu", "beta0", "beta1", "sigma"))

plot(fit_HtWt30_robust, plotfun = "trace",
     pars = c("nu", "beta0", "beta1", "sigma"))

Summarizing the model

print(fit_HtWt30_robust, pars = c("nu", "beta0", "beta1", "sigma"))

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean    sd   2.5%    25%    50%    75%  97.5%
nu     16.17    0.49 11.92   2.95   7.17  12.86  21.76  47.73
beta0 152.19    0.13  4.72 143.17 149.05 152.09 155.14 161.87
beta1   4.31    0.04  1.21   1.85   3.53   4.33   5.12   6.55
sigma  24.05    0.19  4.19  16.65  21.08  23.79  26.67  32.84
      n_eff Rhat
nu      584    1
beta0  1279    1
beta1  1191    1
sigma   481    1

Samples were drawn using NUTS(diag_e) at Fri May 12 12:08:46 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Visualizing the model

pairs(fit_HtWt30_robust,
           pars = c("nu", "beta0", "beta1", "sigma"))

stan_dens(fit_HtWt30_robust,
     pars = c("nu", "beta0", "beta1", "sigma"))

plot(fit_HtWt30_robust, pars = c("nu", "beta0", "beta1", "sigma"))

ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)

The model with 300 observations

fit_HtWt300_robust <- sampling(HtWt_robust,
                               data = stan_data_300,
                               seed = 11111,
                               refresh = 0)

Warning: There were 1 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.Warning: Examine the pairs() plot to diagnose sampling problems

plot(fit_HtWt300_robust, plotfun = "ac",
     pars = c("nu", "beta0", "beta1", "sigma"))

plot(fit_HtWt300_robust, plotfun = "trace",
     pars = c("nu", "beta0", "beta1", "sigma"))

print(fit_HtWt300_robust,
      pars = c("nu", "beta0", "beta1", "sigma"))

Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5%
nu      5.59    0.05 1.79   3.27   4.34   5.25   6.39  10.11
beta0 156.24    0.03 1.67 153.03 155.10 156.22 157.36 159.52
beta1   4.43    0.01 0.41   3.60   4.15   4.43   4.70   5.22
sigma  23.95    0.04 1.66  20.84  22.77  23.92  25.01  27.46
      n_eff Rhat
nu     1142    1
beta0  3146    1
beta1  3097    1
sigma  1656    1

Samples were drawn using NUTS(diag_e) at Fri May 12 12:18:30 2023.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

pairs(fit_HtWt300_robust,
      pars = c("nu", "beta0", "beta1", "sigma"))

plot(fit_HtWt300_robust, plotfun = "dens",
     pars = c("nu", "beta0", "beta1", "sigma"))

plot(fit_HtWt300_robust,
     pars = c("nu", "beta0", "beta1", "sigma"))

ci_level: 0.8 (80% intervals)
outer_level: 0.95 (95% intervals)

Posterior predictive check

We’ll just do this for the example with 30 observations.

Extract samples

samples_HtWt30 <- tidy_draws(fit_HtWt30_robust)
samples_HtWt30

y_sim_post <- samples_HtWt30 %>%
    select(starts_with("y_sim")) %>%
    as.matrix()

Plot results

ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
    geom_point() +
    geom_abline(data = samples_HtWt30,
                aes(intercept = beta0, slope = beta1),
                alpha = 0.01) +
    geom_abline(data = samples_HtWt30,
                aes(intercept = mean(beta0), slope = mean(beta1)),
                color = "blue", size = 2)

Compare to standard linear regression.

ggplot(HtWtData30, aes(y = weight, x = height_mc)) +
    geom_point() +
    geom_abline(data = samples_HtWt30,
                aes(intercept = beta0, slope = beta1),
                alpha = 0.01) +
    geom_abline(data = samples_HtWt30,
                aes(intercept = mean(beta0), slope = mean(beta1)),
                color = "blue", size = 2) +
    geom_smooth(method = lm, color = "red", size = 2)

ppc_intervals(y = HtWtData30$weight,
              x = HtWtData30$height,
              yrep = y_sim_post)

ppc_hist(y, y_sim_post[1:5, ])

ppc_boxplot(y, y_sim[1:5, ])

ppc_dens(y, y_sim[1:5, ])

ppc_dens_overlay(y, y_sim[1:30, ])

ppc_stat_2d(y, y_sim)